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ABSTRACT 

We have extended the ANTARES code to simulate the coupling of pulsation with 
convection in Cepheid-like variables in an increasingly realistic way, in particular in 
multidimcnsions, 2D at this stage. Present days models of radially pulsating stars 
assume radial symmetry and have the pulsation-convection interaction included via 
model equations containing ad hoc closures and moreover parameters whose values 
are barely known. We intend to construct ever more realistic multidimensional models 
of Cepheids. In the present paper, the first of a series, we describe the basic numer- 
ical approach and how it is motivated by physical properties of these objects which 
are sometimes more, sometimes less obvious. - For the construction of appropriate 
models a polar grid co-moving with the mean radial velocity has been introduced to 
optimize radial resolution throughout the different pulsation phases. The grid is radi- 
ally stretched to account for the change of spatial scales due to vertical stratification 
and a new grid refinement scheme is introduced to resolve the upper, hydrogen ioni- 
sation zone where the gradient of temperature is steepest. We demonstrate that the 
simulations are not conservative when the original weighted essentially non-oscillatory 
method implemented in ANTARES is used and derive a new scheme which allows a 
conservative time evolution. The numerical approximation of diffusion follows the same 
principles. Moreover, the radiative transfer solver has been modified to improve the ef- 
ficiency of calculations on parallel computers. We show that with these improvements 
the ANTARES code can be used for realistic simulations of the convection-pulsation 
interaction in Cepheids. We discuss the properties of several numerical models of this 
kind which include the upper 42% of a Cepheid along its radial coordinate and assume 
different opening angles. The models are suitable for an in-depth study of convection 
and pulsation in these objects. 

Key words: stars: variables: Cepheids — hydrodynamics — convection — methods: 
numerical 



1 INTRODUCTION 

The series of papers introduced by this article deals with the 
construction and analysis of time-dependent, nonlinear nu- 
merical models for classical, radially pulsating variable stars 
(Cepheids, RR Lyr,...) in multiple dimensions. The models 
are set up to investigate the pulsation-convection interaction 
and, eventually, nonradial pulsations in such stars. In the 
present paper we describe the extension of the ANTARES 
code ( Muthsam et al.|2010 I necessary for such an investiga- 
tion. We additionally discuss effects of resolution and other 
issues. A subsequent paper ( Muthsam et al.|2012 l, paper II, 
will deal with an investigation of the physical results. 



* E-mail: eva.mundprecht@univie.ac.at 
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The following considerations have motivated this work. 
One of the early successes of numerical modelling of the dy- 
namics of astrophysical objects were simulations of Cepheid 
(or RR Lyr, W Vir,...) pulsations. It was one of the few time- 
dependent nonlinear problems in astrophysics where the ap- 
proximation of perfect spherical symmetry was sufficiently 
good while at the same time observations at the dynamical 
timescale were readily available. Even decades earlier sys- 
tematic trends of the shapes of the lightcurve with period, 
the famous Hertzsprung progression in classical Cepheids, 



had been established (Hertzsprung 19261. Interest in the 



subject continued to be kindled by the importance of these 
variable stars for the understanding of stellar evolution, in 
particular when discrepancies between masses derived from 
stellar evolution considerations and from pulsation proper- 



ties were found (Cox 1980 and Keller 2008). The decisive 
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role of Cepheids in gauging the cosmic distance scale is well 
known (see the extensive discussion in [de Grijs]|201l| e.g.) 
and the precise calibration of relations between mass, lumi- 
nosity, and metallicity for different types of classical vari- 
ables is of general astrophysical interest. 

Linear stability considerations suggested the excitations 
of the pulsations through the ionisation of hydrogen and, 



Again in a review, Marconi J2009I) emphasized the proper 



in particular, helium (Zhevakin 1953 Cox 1960 Baker & 
|Kippenhahn|1962| . The first nonlinear, time-dependent cal- 
culations soon followed. For pioneering nonlinear work see 



Christy 
al.| ( |1966 



19621, Aleshin (19641, Christy (19641, and Cox et 



Spherical symmetry allowed the problem to be treated 
computationally in one spatial dimension and made it ac- 
cessible to the computing resources of that period. The sub- 
sequent decades, essentially till now, have witnessed sub- 
stantial refinements of computational methods and input 
physics. Let us mention the inclusion of radiative transfer for 
the outer zones instead of the diffusion approximation, bet- 
ter numerics including adaptive grids, more realistic opac- 
ities and the like. There now exists a considerable number 
of codes of that kind with various properties such as de- 



scribed in |Fokin| (|1990 


l,|Dorfi & Feuchtinger|(|1991|),|Bono 


& Stellingwerf (19941, 


Buchler et al. ( 1997 1 and Smolec & 


Moskalik (20081. 



It soon became apparent that in order to properly repre- 
sent important observational facts such as the red edge of the 
Cepheid instability strip convection had to be included and 
that attempts to accomplish this by incorporating standard 



mixing length theory into a pulsation code failed (Tuggle 



fc Iben|l973| ). Consequently, practically all groups working 
on the nonlinear hydrodynamics of pulsating stars included 
convection models which from the outset were meant to ap- 
ply to a changing environment, essentially extensions of the 
mixing length approach or models for the time evolution of 
the turbulent kinetic energy (TKE). 

While this improved matters in many respects, it cer- 
tainly did not solve the problem of constructing an accurate, 
predictive model of the convection-pulsation interaction in 
these stars. Since such models cannot be derived from ba- 
sic principles alone, a variety of recipes for modelling the 
convective flux is in use (and is discussed in the next pa- 
per of this series). The same holds for many other physi- 
cal quantities appearing in these models. Many closure con- 
stants, up to eight, "the a's", are needed to fully specify 
the model equations currently used for the one-dimensional 
calculations (see the impressive compilation in |Buchler fc| 
Kollath 2000). Their values are unknown and must be de- 
termined somehow. Common ways to do so include adjust- 
ing the parameters to match some observable quantities or 
simply guessing them. Consequently, even if observed prop- 
erties are reproduced, the extent to which this success can 
be justifiably ascribed to the model and does not constitute 
a mere result of the fitting procedure remains doubtful. It 
most likely lacks the deeper physical correctness intended to 
achieve with unknown consequences for further interpreta- 
tions based on such a model. 

The state of knowledge has recently been discussed in 



the review by Buchler (20091. There, he addresses the use 
of time dependent mixing length theory as the "dominant 



inclusion of convection for modelling RR Lyr stars. 

Recently, activities have commenced to address this is- 



Gastine & Dintrans (20111 thoroughly analyse phys- 



ically simplified 2D models to elucidate basic mechanisms 
of the pulsation-convection interaction (see also literature 



cited therein). Geroux & Deupreel (20111 have reassumed 



very early work by Deupree ( 19801). Their calculations are 



in 3D, albeit with quite low resolution at present. Various 
physical simplifications have been made in this work (no re- 
alistic microphysics, adiabatic stratification), but more real- 
istic microphysics, opacities, and an energy equation based 
on the diffusion approximation are intended to be included 
in future investigations. 

This is the context in which we present the extension of 
the ANTARES code for the type of research at hand. We ex- 
tend the main features of the code (radiation-hydrodynamics 
with radiative transfer equations in the outer regions, real- 
istic microphysics and opacities, high-resolution numerics, 
optional grid-refinement) as required for work on nonlinear 
stellar pulsation. We furthermore add a moving grid (mov- 
ing at the upper boundary with the mean vertical surface 
velocity) . While the main part of ANTARES is designed for 
the ID, 2D, and 3D case, for our present work only the ID 
and 2D case have been implemented. This is due to the fact 
that the extremely steep gradients near the hydrogen ionisa- 
tion zone make even an adequately resolved 2D calculation 
a computationally very expensive enterprise. Undoubtedly, 
however, in the future we will also include the option of 
performing numerical simulations of this kind in 3D. 

At first sight it may look astonishing that we are go- 
ing to undertake 2D calculations for the time being. After 
all, meaningful 3D models of solar granulation exist since 



decades (Nordlund 1982). For Cepheids, however, the cir- 
cumstances are different. We find a few remarks on the 2D 



model atmospheres of Cepheids in Freytag et al. ( 2012 I. The 
difficulties are basically the same as those encountered in 
the simulation of A-type stars and extensively discussed 



for them by Kupka, Ballot, & Muthsam (20091 and, for 
Cepheids, in Sect. [3j of the present article where also a few 
remarks on computational demands are provided. - The 3D 
nature of the simulations presented in Geroux fc Deupree| 
(20111 does not contradict those remarks, since their work 



considers an adiabatic test case whose idealised microphysics 
lacks the properties which ultimately lead to the aforemen- 
tioned problems. 

With the goal of modelling such stars in multidimen- 
sions with an ever increasing degree of realism we have ex- 
tended the ANTARES code and continue to do so, maintain- 
ing its basic design principles and ingredients as described 



Muthsam et al. (20101, namely 



deficiency" in Cepheid modelling (see also Buchler (1997 1 ) 



• high resolution numerics of the ENO variety 

• realistic equation of state and microphysics 

• radiative transfer in the diffusion approximation or us- 
ing the static radiative transfer equations 

• parallelization based on MPI and OpenMP. 

For our present specific purpose extensions are neces- 
sary which are either obvious from the outset or which ex- 
perience has taught to be indispensable, namely 

• the use of polar coordinates 
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• a moving grid to follow the pulsations of the star 

• a radially stretched grid for coping with the vast dif- 
ferences in length scales between the stellar surface and the 
interior 

• adaption of the treatment of hydrodynamics and radia- 
tive transfer to these changed paradigms 

• specific forms of grid refinement to properly cope with 
the fast spatial variation of quantities in very localized re- 
gions 

• considerations of how to formulate the numerical 
scheme given the requirements of the objects envisaged for 
research. 

The aim of the paper is, therefore, 

• to discuss just those numerical considerations and how 
they are brought about by the physics of pulsating stars 

• and to provide an overview of what is accessible to mul- 
tidimensional numerical modelling of Cepheids and related 
variables presently or in the foreseeable future. 

For that purpose we have organised the remaining sec- 
tions as follows. We first present the equations of radiative 
hydrodynamics on a stretched, co-moving and polar grid 
and some numerical methods adapted to this setting. In 
Sect. [3] we provide the parameters of our model and inves- 
tigate the restrictions resulting from the steep gradient in 
the H-ionisation zone, while in Sect. [4] we analyse the ef- 
fects of resolution. Finally, we discuss the usefulness of the 
calculated simulations, their strength and weaknesses, and 
possible improvements in Sect. [5] before we provide our con- 
clusions in Sect. [6] 



9e =-V- 
dt 



- (e + p) - eu g 
P 



where the term containing 



I 

— ■ a 
P 



2p + 



+ g-i + < 



-Irlip+Iiple cot f 



y — p cot 9 



(3) 



(4) 



originates from introducing spherical coordinates, a is the 
viscous stress tensor, p denotes density, I = (I r , I v ,Ie) is the 
momentum density, e the total energy density, p the pres- 
sure, and g the gravitational acceleration. The grid velocity 
at the boundaries is set to the horizontal fluid velocity aver- 
age at the top, u g \top = Wr|top, and zero at the bottom. For 
the intermediate points u g (i) = (n - r bot ) / (r top - r bot ). 
The radiative heating rate Q ra( j is computed by solving the 
radiative transfer equation (RTE) in optically thin regions, 
and by the diffusion approximation —V • (xVT) elsewhere. 
X is the radiative conductivity. The system is closed by an 
equation of state, which describes the relation between the 
thermodynamical quantities and depends on the physical 
properties of the fluid. Since the simple thermodynamical 
relations for an ideal gas are not applicable to a partially 
ionised gas, realistic microphysics is included by the LLNL 



equation of state (Rogers et al. 1996 Rogers & Nayfonov 



2002[ ) an d OPAL opacities ([Iglesias fc Rogers|1996| |. For the 



RTE the |Ferguson et al.| ( |2005[ ) low-temperature Rosseland 
opacities for grey radiative transfer are used to extend the 
accessible temperature range. 



2 NUMERICAL METHODS 



2.2 The grid structure 



The simulations of radially pulsating stars are performed on 
a stretched, polar, and moving grid. We first introduce the 
basic conservation laws governing the time evolution of our 
numerical models in Sect. 



2.1 



2.2 



and define subsequently the 
. In Sect. 12.31 the necessary 



computational grid (Sect, 
alterations for the radiative transfer equation (RTE) are ex- 
plained in detail. They are not only required for realistic sim- 
ulations of layers near the stellar surface, but also alleviate 
the restrictions on time stepping (see Sect. 3.2 1. In Sect. 2.4 a 



new set of ENO-coefficients is presented that ensures conser- 
vation of the hydrodynamic quantities on a stretched polar 
grid. Similar adaptations for discretising diffusive fluxes are 
discussed in Sect. [275] while Sect. [2~6] explains our optional 
subgrid scale modelling. Finally, Sect. |2. 7] explains the setup 
of initial and boundary conditions of the simulations. 

2.1 Conservation equations 

Including the grid velocity in the conservation equations 
for mass, momentum, and total (kinetic plus thermal) en- 
ergy we obtain as our set of dynamical equations 



dp 
dt 



= -V ■ [/- put. 



97 
dt 



= -V 



n 

p 



Iu a 



a + p ■ Id 



source 
H h pg, 



(1) 



(2) 



The computational domain is equipped with a 3D spherical 
grid, which is used for 2D simulations and for ID simula- 
tions, such that one cell covers one entire shell at a given 
radius. In radial direction N r grid points plus 4 ghost cells 
at each boundary are used. The values of N r etc. for the dif- 
ferent models can be found in Table [I] The r-range covers 
r G [rbot, rtop], where rbot is fixed and r top varies with time. 
The grid is stretched in radial direction by a factor q. The 
mesh sizes are Ar i+1 —qAri varying from top to bottom. 
Thus, the numerical grid becomes 



: r bo t + 



(r top — r bo t) • 



q"<* — 1 

Cell boundary values are defined as r i+ i = n 



(•5) 



and r- 1 = ti + %/ f Ar '_ 1 . In azimuthal direction the angle 
2 t+v*? 

ip covers tfj±i G [— s jjr t ,+ ie jr 
i ^ et . The numerical grid is given by 



and the mesh size is Aip = 



^tOt . / ■ 1 \ A 

<Pi = + U - 2> A V- 



(6) 



Cell boundary values are defined as <Pj±i = <fj ± -j 2 . The 
physical distance between two adjacent points in azimuthal 
direction is computed as Ayit — r-iAip sin Ok- The polar (co- 
latitudinal) angle covers the range 6 k± i G [ n ~ 2 tot > 7r+ 2 t ° t ] > 
the mesh size is AO — ^f^- The numerical grid is given by 

- tot - AO 



k = kAO 



© xxxx RAS, MNRAS 000,[T]{T5] 



4 E. Mundprecht, H.J. Muthsam and F. Kupka 



model 


aperture 


grid points 


stretching 


subgrid 


grid 


radial 




refined radial 


nr. 


angle 


radial 


polar 


factor 


modelling 


refinement 


cell size 




cell size 


1 


10° 


510 


800 


1.011 


no 


no 


0.47 . . . 


124 Mm 




2 


3° 


510 


300 


1.011 


yes 


no 


0.47 . . . 


124 Mm 




3 


3° 


510 


300 


1.011 


yes 


yes 


0.47 . . . 


124 Mm 


0.32 ...0.80 Mm 


i 


1° 


800 


300 


1.007 


yes 


no 


0.29 . . . 


79 Mm 





Table 1. Grid parameters for the different numerical Cepheid models discussed in this paper. 



Table 2. Quadrature formula Carlson A4 for integration over 
a sphere: ray directions in the first octant. 



Cell boundary values are defined as 9 k± i = 9k ± — . The 

physical distance between two adjacent points in colatitudi- 

nal direction is computed as Az 4 = r^Af?. 

For 2D simulations there is just one cell of size ^tot = 

AO — 7r in colatitudinal direction and the grid reduces to 

6i = 0, 6*i = ~, and 63 = tv and for ID simulations in 

2 A 2 

addition the azimuthal grid becomes ip 1 = — tt, ipi = 0, and 

ip3 = 7T with <£> t ot = A</3 = 2-7T, so that one cell represents one 

shell at a given radius. In 2D simulations we hence consider 

sectors located in the equatorial plane of the sphere as our 

simulation boxes while in ID the simulation boxes are radial 

columns. In ANTARES the :r-direction of the Cartesian grid 

points inward, thus for compatibility with various routines 

Xi = ftop — Ti is used and for visualization purpose there 

is a Cartesian output grid. There, the x-coordinates of the 

radius at p> — and = 5 appear negative which has to be 

considered when evaluating and visualizing the output data. 



2.3 The radiative transfer equation (RTE) 

To perform a realistic simulation of layers near the stellar 
surface the nontrivial exchange of energy between gas and 
radiation is included via the RTE, while the diffusion ap- 
proximation is used for the deeper, optically thick layers to 
reduce computing time. The RTE in ID 



di 

^JT = 1 

OT 



S 



(8) 



is solved along single rays via the short characteristics 



method of Kunasz & Auer ( 1988 1. Here, I denotes the inten- 



sity, S the source function (taken to be the Planck function), 
and t the vertical optical depth scale. Finally, fi — cos 6, 
where 9 is the polar angle. In 2D either 12 or 24 ray direc- 
tions are chosen according to the angular quadrature formu- 



lae of types A4 or A6 of Carlson ( 1963 1 , and the directions in 



each quadrant are arranged in a triangular pattern (Tables 
[2]and|3|. This is a reduction of a factor of two compared 
to the 3D case, a result which follows from the additional 
symmetry. 

For each ray the points of entrance and exit plus the 



1 




7i 


z» 


1 




pT 

V 15 


fl3 
V 15 


2 


nr 

V 15 


pf 

V 15 


pf 
V 15 


3 


fT~ 

V 15 


fl3 
V 15 


fT~ 

V 15 


4 


pf 

V 15 


V 15 


AT" 

V 15 


5 


pf 
V 15 


rv 

V 15 


PC 

V 15 


6 


/l3 
V 15 


pT 

V 15 





Table 3. Quadrature formula Carlson A6 for integration over 
a sphere: ray directions in the first octant. 
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Figure 1. One step in the iterative solution of the RTE. The 
rays enter at U. There, at the center P, and at the exit point Q the 
opacity has to be determined to perform the numerical integration 
on the right hand side of Eq. |9]| and to get the intensity / at P 
for the next step. 



corresponding distances are determined (Fig. [T]). Since the 
grid moves this has to be redone at every time step. The one 
dimensional RTE is solved along each ray. 

To determine the incoming intensity along this ray di- 
rection one has to interpolate the values of the required phys- 
ical quantities to the points U and Q from the neighbouring 
grid points. This is achieved by linear interpolation. Conse- 
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k 


r 


j=o 


j=l 


j=2 


j=3 


j=4 


1 


-1 


1 















1 










2 


-1 


3/2 


-1/2 













1/2 


1/2 










1 


-1/2 


3/2 








3 


-1 


15/8 


-5/4 


3/8 











3/8 


3/4 


-1/8 








1 


-1/8 


3/4 


3/8 








2 


3/8 


-5/4 


15/8 






4 


-1 


35/16 


-35/16 


21/16 


-5/16 









5/16 


15/16 


-5/16 


1/16 






1 


-1/16 


9/16 


9/16 


-1/16 






2 


1/16 


-5/16 


15/16 


5/16 






3 


-5/16 


21/16 


-35/16 


35/16 




5 


-1 


315/128 


-105/32 


189/64 


-45/32 


35/128 







35/128 


35/32 


-35/64 


7/32 


-5/128 




1 


-5/128 


15/32 


45/64 


-5/32 


3/128 




2 


3/128 


-5/32 


45/64 


15/32 


-5/128 




3 


-5/128 


7/32 


-35/64 


35/32 


35/128 




4 


35/128 


-45/32 


189/64 


-105/32 


315/128 



, 5, j = 0, 



, 4 on 



an equidistant grid. 



quently, in two dimensions 4 points are used for all quantities 
except for the intensity where 2 points are sufficient. Then 
we evaluate the equation 

r T p 

I (t p ) = I (t v ) exp (r„ - t p ) + S (t) exp (t - r p ) dt (9) 

Jtu 

numerically to get 7 (tp) from a quadrature rule proposed 



by Olson & Kunasz (19871. This procedure is repeated re- 
cursively, because after the first step one obtains just the 
intensity on a single new point. 

The mean intensity J is computed by solid angle inte- 
gration. For A4 the weights are identical, 



E 



(10) 



Finally, the radiative heating rate for grey radiative transfer 
is given by 



Qrad = 471-pX (J - S) , 



(11) 



as in the case of a non-moving grid (details on how this is 
handled within ANTARES can be found in Muthsam et al.l 



(20101). 



2.4 WENO schemes on stretched and co-moving 
polar grids 

The differential equations solved by the ANTARES code can 
also be written as a system of conservation laws of the form 



dtU = — V • F (U) + S. 



(12) 



Here, F is the vector-valued flux function and S designates 
(physical or geometrical) source terms. A discretisation of 
the physical domain results in a centre grid n and a bound- 
ary grid r i+ i, where the function Fi(U) is computed at 



the cell centre and interpolated to the cell boundaries. If a 
simple polynomial interpolation from the centre grid to the 
boundary grid is used, shocks or discontinuities are either 
smoothed out or begin to oscillate. To avoid these problems 



various essentially non-oscillatory (ENO) schemes (Harten 
etaL] jT987| , |Shu fc Osher| ( p88| , piu| (m7\ , [Fedkiw et al. 
1 1998| )) have been implemented into the ANTARES code, 
where interpolation is done by adaptive upwinding stencils. 
In the ANTARES code the spatial discretisation is done for 
each direction separately. However, note that for a polar grid 
the radial and angular direction are different. In radial di- 
rection the grid is stretched while it is equidistant in angular 
direction. The basic idea of characteristic numerical schemes 
is to transform this nonlinear system into a system of nearly 
independent scalar equations, a typical one of the form 



«t + f(u)x = 



(13) 



(possibly with an inhomogeneity at the right hand side) 
and discretise each scalar equation independently and then 
transform the discretised system back into the original vari- 
ables. Thus, for hyperbolic systems the set of conservation 
laws is first transformed to the eigensystem, where the equa- 
tions decouple and the upwinding directions can be chosen. 
For the moving grid the eigenvectors are left unchanged as 
the grid velocity u g only enters into the eigenvalues: 



(14) 



Schemes of the ENO variety as first introduced by 



/ 


u x 


- Ug 


^ \ 




u x 


- U g 






U x 


- Ug 




V 


U x 


- Ug 


+ c ) 



Harten, Engquist and Chakravarthy (Harten et al. (1987 1 ) 
are based on cell averages as follows. Given the cell aver- 
ages /j = / (n) of a function we want to find a numerical 
flux function f i+ i = f (fi-r, ■ ■ ■ , fi+s) such that the flux 
difference approximates the derivative /' (r;) where / is suf- 
ficiently smooth. The mapping from the given cell averages 
{fjj in the stencil S (i) to the values f~ +1 / 2 1& unear - There- 
fore, there exist constants c r j and c r j depending on the left 
shift r of the stencil S (i) such that 



fc-i 

fi + l/2 ~ C rj.fi-r + j 

j=0 



and 



-1/2 



fc-1 
J=0 



-r+j 



(15) 



(16) 



where c r j = c r _i j . Near discontinuities in the solution of hy- 
perbolic conservation laws oscillations can occur because the 
stencil contains those discontinuities. Therefore, an adaptive 
stencil S (i) is chosen for the interpolation of the cell bound- 
ary fluxes where the left shift r changes with the location 
r%. After all the main idea of the ENO approximation is to 
exclude cells containing discontinuities from the stencil S (i). 



When using WENO (Liu et al. (19941, Jiang fc Shu 



( |1996 1) instead of choosing a single stencil for the interpola- 
tion polynomial in the ENO reconstruction, a convex com- 
bination of all candidates is used to achieve the essentially 
non-oscillatory property. A weighted combination of adap- 
tive stencils gives high order polynomial approximations of 



© xxxx RAS, MNRAS 000,[T]{T5] 



6 E. Mundprecht, H.J. Muthsam and F. Kupka 



k 


r 


j=o 


j = l 


j=2 


3 


-1 


a 6 +2a 5 +3a 4 + 3a 3 +3a 2 +2a + l 


a 4 + a 3 +a 2 +a+l 


a 2 +a+l 


a 2 (a 4 +2a 3 +2a 2 +2a + l) 


a 4 (a 2 +2a + l) 


a 4 (a 4 + 2a 3 +2a 2 +2a+l) 







a 2 (a 2 +a + l) 


a 2 +a+l 


1 




a 4 +2o 3 +2a 2 +2a+l 


a(a 2 +2a+l) 


a(a 4 +2a 3 +2a 2 +2a + l) 




1 


a 5 


a(a 2 + a + l) 


a 2 + a+l 




a + 2a 3 + 2a' 2 +2a + l 


a+2a + l 


a 4 +2a a +2a^+2a + l 




2 


a 6 (a 2 +a + l) 


a 2 (a 4 + a 3 + a 2 +a + l) 


a 6 +2a 5 +3a 4 + 3a 3 +3a 2 +2a+l 




a 4 +2a 3 +2a 2 +2a+l 


a 2 +2a + l 


a 4 +2a 3 +2a 2 +2a + l 



Table 5. The constants C r j for k = 3, j = 0, . . . , 2 on the stretched grid. For better readability, a = s fq. 



k=5, r=2 



a 16 (a 2 +a + l) 

a 14 + 3a 13 +5a 12 +8a 11 + lla 1 ° + 13a 9 + 15a 8 + 16a 7 + 15a 6 + 13a s > + lla 4 +8a 3 +5a 2 +3ti + l 

a 8 (a 4 + a 3 +a 2 +a + l) 
_ a 8 +3a 7 +4a f >+5a 5 +6a 4 + 5a :i +4^+3(1+1 
a 2 (a 8 +3a 7 +6a 6 +8a 5 +9a 4 + 8a 3 +6a 2 +3a+l) 
a s +4a 7 + 8a 8 + 12a 5 + 14a 4 + 12a :i +8a :4 +4a + l 

a 6 +2n B +3a 4 + 3a 3 +3a 2 +2a + l 

a(a 8 +3a 7 +4a 6 + 5a B +6a 4 + 5a 3 +4a 2 +3a+l) 

a 4 +a 3 +a 2 +a + l 

" a(a 14 + 3a 13 +5a 12 +8a 11 +lla 10 + 13a 9 + 15a 8 + 16a 7 + 15a 6 + 13a 5 + lla 4 +8a 3 +5ci 2 +3a + l) 



Table 6. The constants C^j for k = 5 and r = 2. 



the divergence term V • F (U) in smooth regions. Across dis- 
continuities the smoothest stencil is chosen even though it 
is of lower order. Instead of performing a 2k — 1 order ENO 
scheme using stencils of that order a combination of k sten- 
cils of order k is used to obtain the final accuracy 2k— 1. For 
ANTARES the orders 3 and 5 have been considered. For a 
fe-th order ENO scheme there are k candidate stencils 

S r (i) = {Xi-r, ■ ■ ■ ,Xi-r+k-l} , r = 0, . . . , k - 1 (17) 

which produce k different reconstructions of the value f i+ i ■ 

k-l 



f i + 1/2 — E Cr lfi-r+j> r — 0, 
3=0 



,k 



(18) 



r I 

The convex combination of the values / i+1 /2 ^ or l ' ne WENO 
approach is 



/i+i/2 = J2 Wr ^ 



(r) 

+1/2' 



r = 0, 1 (19) 

and is used as a new approximation for /i+1/2 • For the 
weights uv for all r and X^r=o "V = 1 must be true. 
For a smooth function / (x) there are constants d r so that 



/i+1/2 = E rf ^+i/ 2 = / + ° ( A ^ fc_1 ) 

and constants d r so that Xlito dr = 1 and 



(20) 



/i-l/2 — E drf i+l 



/2 



1/2. 



+ O Ar 



(21) 



In Liu et al. (19941 weights of the form 



Ek-1 
s=0 ' 



(22) 



are proposed, e > is the machine ac- 



with a r = - {E+0r) , 
curacy, which is introduced here to ensure that the denom- 
inator does not become zero. /3 r are smoothness indicators 
of the stencil supposed to be zero when a discontinuity is 
contained in the stencil. A robust choice of smoothness in- 
dicators is defined by 



i+1/2 



E 



A 2 

Ar,- 



8'pr (x) 

dx l 



dx 



(23) 



where p r (x) is the reconstruction polynomial on the stencil 
S r (i)- The smoothness indicators f3 r are a measure for the 
total variation in the interval Ii. 

From a different point of view another set of approxi- 
mation coefficients C r j (cf. Tables [4]j6| , of weights w r (cf. 
Table [7]), of their approximations D r and D r and of smooth- 
ness indicators, denoted by being capitalised, can be derived. 
In this case the smoothness indicators for k — 3 are: 

A) = Vo 2 (/, -(q + 1) fi+i + fi+2) 2 + 

(Ziofi + -£20/1+1 + Zzofi+2) 2 

0i= if (7i-(g + i)/i+i + /i+2) 2 + 

(■£11/1-1 + Z21.fi + -£31/1+1) 
02 = Y? (/, -(q + 1) /i+i + / i+2 ) 2 + 
(■£12/1-2 + -£22/1-1 + -£32/1) 



(24) 



where F and Zi are given in Table [8] These are also available 
in the ANTARES code in addition to those from ( 23 1 which 
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D 

D-2 



Do 
Dl 





^a s +2a 7 +2a (i +2a f > + 2a 4 + 2a a - 


r a 2 +a+l 




a 3 (a 4 +a 3 +a 2 +a + l) 






a s +a 7 +a 5 +2a 4 +a a +a + l 






a 10 






-a 8 +2a 7 +2a 6 +2a 5 + 2a 4 + 2a 3 - 


-a 2 +a+l 


1 




-a 8 +2a 7 +2a 6 +2a 5 + 2a 4 + 2a 3 - 


-a 2 +a+l 




a(a 4 +a 3 +a 2 +a + l) 






a 8 +a 7 +a 5 +2a 4 +a a +a + l 




„ll) _!_„!)_ 


a 6 (a 4 + a 3 +a 2 +a+l) 





Table 7. Weights D r for f~ +1/2 and D r for /+_ 1/2 for fc = 3 



Z 2 j 



Z 3i 






\/l3 


a 2 +a+l 


a 4 + a 3 +a 2 +a+l 


1 


V3a 4 (a 2 +l) 


a 2 (a + l) 


a 4 (a + l) 


a 4 (a+l) 


1 

2 




a 2 


-a 3 + l 


1 


V5(a.2 + 1) 
yi3a 4 


a+1 
a 5 


a(a+l) 
a(a 4 +a 2 +a + l) 


a(a+l) 
a(a 2 +a + l) 


^3(a 2 + l) 


a + 1 


a+1 


a+1 



Table 8. Coefficients for evaluating f}j, j = 0, 1,2, for k 
Equation [24] 



are based on the approach ( 19 H 22 l and are required for 



stable time integrations of our Cepheid models, as we show 
below. In the second approach instead of viewing f (ri) as 
cell averages and interpolating the primitive function, the 
flux function / (ji) itself is interpolated to the boundaries, 
leading to a different set of coefficients (Table j4j). In our case 
a combination of stencils of length 3 with coefficients C r j 
from Table [5] yields the stencil of length 5 with coefficients 
Cij in Tableplin smooth regions. This approach was chosen 
because we use a polar grid and although the spacing is 
that of a stretched grid in radial direction, the cell-volumes 
depend also on the radius as does the shape of the cells. 
The approach via exact cell volumes ( 19 H 23 1 is not 



advisable in our case since we want to ensure conservation 
of the basic variables. Though the difference in area size be- 
tween a rectangular cell and a polar cell is very small the 
difference over the whole grid is not and the ensuing error 
cannot be neglected. When using the original coefficients in 
a one-dimensional simulation e.g. energy is not preserved 
over a prolonged period of time: during the first 48 days 
of a Cepheid simulation it increases or decreases depending 
on various conditions. (The parameters for this Cepheid are 
provided Sect. 



Fig.[2f, 



3.11 



One run shows increasing values (see 
but for a slightly larger domain and with artificial 
diffusivities included we obtain decreasing values. With the 
new set of coefficients all the conserved variables are stable 
(Fig. [3} when averaged over several pulsation cycles. This 
comparison also demonstrates that a non-conservative dis- 
cretisation method is unacceptable for a realistic numerical 
simulation of Cepheids. 

Note that the formally fifth order of the WEN05 



scheme ( Shu 2003 1 does no longer hold with the new ap- 
proach. However, lack of energy conservation is not an ac- 
ceptable alternative for locally higher order in smooth re- 
gions. Around discontinuities both old and new coefficient 
sets yield methods of second order which is also the spa- 




Figure 2. Cepheid simulation with interpolation of the prim- 
itive function with cell averages in all directions: time evolution 
of radius (red) and total energy (green). Abscissa: time in days, 
left ordinate: radius on a linear scale relative to its initial value, 
right ordinate: total energy within the sphere in erg. 




Figure 3. Cepheid simulation with interpolation of the flux 
function and the new scheme used in all directions: time evolution 
of radius (red) and total energy (green). Abscissa: time in days, 
left ordinate: radius on a linear scale relative to its initial value, 
right ordinate: total energy within the sphere in erg. 



tial order of accuracy achieved with the new set in smooth 
regions. Moreover, for the one-dimensional linear advection 
equation with constant velocity the error coefficient is one 
quarter of the size of standard second order methods such 
as the leap-frog scheme. Leap-frog time integration, how- 
ever, can cope with shocks only after modifications that de- 
grade their accuracy even further. The main computational 
expense of the WENO approach is the calculation of the 
fluxes at each cell boundary. This occurs only once, whence 
the overhead due to larger interpolation stencils is negligi- 
ble, as they merely use already computed information. The 
only important restriction originating from wider stencils 
in comparison with traditional second order shock captur- 
ing methods are somewhat larger minimum domain sizes for 
the domain decomposition approach used in the MPI paral- 
lelization in ANTARES. 
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_ 3 + 6a+<)a z + 10ri +lla +l()a ,J +l)a °+(ki' + ia*+2a J + a L 



a 2 (l+a) 2 (l+2a 2 +2a 4 + a 6 ) 

3 + 4a+5a 2 +6a 3 + 7a 4 + 6a 5 +6a 6 +4a 7 +3a 8 +2a 9 +a 1 
a 6( 1 + a )2( 1+a 2) 

3+4a + 5a 2 +4a 3 +4a 4 + 4a 5 +4a 6 + 2a 7 + a 8 
a 8 (l + a) 2 (l + a 2 ) 

3+4a + 5a 2 +4a 3 +4a 4 + 2a 5 +a 6 
a 8 (l+a) 2 (l + 2o 2 +2a 4 + a 6) 



71 

" 24 



47 

8 



23 
24 



3a a +4a a +5a'' + 4a'+4a g + 2a J +< 
(l + a) 2 (l + 2a 2 +2a 4 + a 6) 

-2a-a 2 +2a 4 +2a +3a 6 +2a 7 + a 
a 2 (l + a) 2 (l + a 2 ) 

2a-a 6 
a 4 (l + a) 2 (l + a 2 ) 

2a-a 4 

a 4 (l + a) 2 (l + 2a 2 +2a 4 +a<5) 



23 
24 



1 

" 24 



(l+n) 2 (l+2a 2 +2a 4 +a6) 
(2 a 3 -a 4 )(l + a + a 2 ) 2 

(l+ a ) 2 (l+a 2 ) 
(-l + 2a)(l + a + a 2 ) 2 
(l + a) 2 (l + a 2 ) 

(-l + 2a 3 ) 

(l + a) 2 (l + 2a 2 +2a 4 + a 6) 



J_ 
24 



J_ 

' 24 



-2 


-a 12 +2a 15 
(l + a) 2 (l+2a 2 +2a 4 +a6) 


1 

24 


-1 


-a 6 +2a 11 
(l + a) 2 (l + a 2 ) 


1 

8 





a 2 +2a 3 +3a 4 + 2a 5 + 2a 6 -a 8 -2a 9 


7 


(l + a) 2 (l + a 2 ) 


8 


1 


a 2 +2a 3 +4a 4 +4a 5 +5a 6 +4a 7 + 3a 8 


23 


(l + a) 2 (l + 2a 2 +2a 4 +a<i) 


24 



a 1 "+2a 1 -'+4a 1D +4a ll +5a^+4a 1 



(l + a) 2 (l + 2a 2 +2a 4 + a | 3) 

a 8 +2a 9 +4a 10 +4a 11 +4a 12 +4a 13 +5a 14 + 4a 15 +3a 16 
(l+a) 2 (l+a 2 ) 

f3a 6 +6a 8 +7a 1Q +5a 12 +3a 14 +2a 5 +4a 7 +6a 9 +6a 11 +4a 13 
(l + a) 2 (l + a 2 ) 

-, a 4 + 2a o +4a 6 + 6a 7 +9a 8 + 10a 9 + lla 10 + 10a 11 +9a 12 +6a 13 +3a 14 







(l + a) 2 (l + 2a 2 +2a 4 + a f 5) 



23 
24 



71 
24 



Table 9. Discretisation of the derivative terms representing dif- 
fusion. Weights for the stencil (i — 2, i — 1, i, i + 1) for the deriva- 
tives at, from top to bottom, i — 2.5, i — 1.5, i — 0.5, i + 0.5, and 
i + 1.5. The derivatives at % — 2.5, i — 1.5, i + 0.5, and i + 1.5 are 
only needed near the boundaries, in the inner regions the weights 
for i — 0.5 are used (cf. Eq. | |26| l). For the equidistant grid these 
terms reduce to the values in the rightmost column. For better 
readability, a = ^fq. 



2.5 Discretisation scheme for diffusion 

To avoid problems with conservation of mass, momentum, 
and energy as exemplified through Fig. [2] the interpolation 
of the diffusive fluxes has to be done in the same manner. 
The new WENO coefficients are given b y Tabl es [5]{8l |HaFl 
|penhofer et al.| ( |2011[ ) and |Koch et aT] ( |2010| ) discussed a 
differencing scheme for diffusive fluxes based on interpolat- 
ing cell averages. For constant diffusivities and equidistantly 
spaced grid points in one dimension in a Cartesian coordi- 
nate system it recovers the well-known approximation 



Ui- 



16 Ui-i +mUi- 16 Ui+i + U t+ 2 



which is of fourth order. However, this cannot be carried over 
to the case of a stretched, co-moving, polar grid for the same 
reasons as already discussed for the WENO scheme. Plain 
staggered mesh interpolation of fluxes can be used instead of 
interpolating the cell averages in which case the derivative 
at i + a — \ is computed as 



S-aUi-i + S-tUi-! + S Ui + SiUi 



for a - 



and 



U, 



n+i — ri 

2. For the equidistant grid this leads to 
24 h 



(25) 



(26) 



Ui 



-27U t + 27U l+1 ~U t+ 2 



24 h 



and the second derivative is then given by 

Ut-2 - 28 Ui-i + 54 Ui - 28 U i+1 + U l+2 



(-Ui) rr = 



24 ft 2 



(27) 



(28) 



(the latter was also discussed in |Koch et al.| ([2010) for Carte- 
sian grids). This is a second order approximation for (Ui) rr 
with an error coefficient half the size of that one of the cen- 
tral second order difference on an equidistant grid in one 
dimension for the heat equation with constant diffusivity. 
As with the WENO scheme, the expensive part in evaluat- 
ing these stencils is the computation of the cell boundary 
fluxes, whence the extra accuracy comes at negligible costs. 
Moreover, the stencils are still smaller than their WENO 
counterparts and thus introduce only small overheads when 
used in ANTARES in parallel mode. This approach can be 
generalised to the stretched and co-moving polar grid in a 
conservative manner. Details are given in Table [9] which dis- 
plays these new stencils for both the interior of the domain 
(third panel) and for points near the (vertical) boundaries. 

2.6 Subgrid scale modelling 

In some of our simulations the grid refinement was aided 
or substituted by subgrid scale modelling. The approach we 
use for this purpose is based on the work of |Smagorinsky| 
(19631 and Lilly (19621. The subgrid-scale stress is written 



3 5ijTkh 



-2pK m Di 



(29) 



where the superscript "a" denotes the anisotropic part of a 
tensor and Da is the strain rate tensor, 



di 



di 



dxj dxi 



1 



dilk 



The eddy viscosity K m is given by 
K m = CA 2 \d\ 



(30) 



(31) 



12 h? 



where A is the filter width, approximated in 2D as 
(Ari Aj/i) 1//2 , and as (AnAyiAzi) 1 ^ 3 in three dimensions, 
and C is a dimensionless coefficient. This coefficient is a 
constant parameter which is commonly expressed through 
the Smagorinsky coefficient c s = C 1 ^ 2 . For our simulations 
we have chosen c s = 0.2 which is the standard value found 
in the literature for 2D simulations (for 3D simulations a 
smaller value of about 0.1 is usually preferred). 
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2.7 Initial and boundary conditions 

We start with a one-dimensional model. For the simulations 
presented here such a model was kindly provided by Giinter 
Houdek. There is no turbulent pressure included, thus the 
starting model set up in ANTARES is purely radiative. To 
accelerate the onset of pulsations the gravitational accelera- 
tion is reduced by 1% during the first few time steps. After 
the pulsation becomes stable and remains so for a consid- 
erable time the one-dimensional model is converted into a 
two-dimensional one by copying the one-dimensional model 
repeatedly in azimuthal direction. Slight random perturba- 
tions are applied to the new lateral momentum variable to 
start the truly 2D flow. 

For closed boundary conditions the density is set to hy- 
drostatic equilibrium at the top. Moreover, the conditions 



| top 

du v . 

~Qf I top 

dT 

w ]top 

are applied at the top, and 



u g | top 






Wr|bot = 



(32) 



(33) 



at the bottom. Since the star is purely radiative at the bot- 
tom the incoming energy transport is adjusted by keeping 
the radial component of radiative flux density kVT at its 
initial value. After applying the boundary conditions, the 
grid is updated by moving the grid coordinates with grid 
velocity Ug. The new grid velocity u g l+ is determined as the 
horizontal average of the vertical velocity at the top. 



3 COMPUTATIONAL CHALLENGES 

In this section we discuss, after presenting the parameters of 
our model, various restrictions on the timestep and their in- 
teraction with grid refinement. For model 1 the timestep 
varied from 0.12 to 0.16 seconds. This leads to approxi- 
mately 2.5 million Runge-Kutta steps for the computation 
of one pulsation cycle. The simulation was carried out in 
MPI-mode on 256 nodes. At continuous operation this al- 
lows the calculation of one pulsation cycle every five days. 
For the better resolved models the computational demands 
were even higher. 



3.1 The model parameters for the Cepheid 
simulations 

For our Cepheid models we assume an effective tempera- 
ture of T c ff = 5125 K, a surface gravity of 92.48 cms -2 (i.e. 
log(g) ~ 1-97), a luminosity L ~ 913 Lq, and a mass of 
5 Mq. For their composition we assume a subsolar metal- 
licity, i.e. a hydrogen mass fraction of X = 0.7 and a metal 
mass fraction of Z = 0.01 (with the composition of Z taken 
from Grevesse & Noels||1993 1. These values are well within 



the parameter range for which Cepheids can be observed 



(see Bono et al. 19941. Through a series of tests we found 



that these model parameters also have the advantage of re- 
ducing the numerical demands in comparison with higher 
luminosities and lower metallicities, for example. 



The stellar radius in our models hence is 26.8 Gm (i.e. 
R ~ 38.5 Rq), of which the outer 11.3 Gm were modelled. 
Three different aperture angles were used: 1°, 3°, and 10°. 
The 3°-models calculations were also performed with a grid 
refinement in the hydrogen-ionisation zone with refinement 
factors of 3 in radial direction and 4 in polar direction. Fur- 
ther details such as the radial stretching factor between ver- 
tically adjacent grid cells (innermost cells are the longest 
ones), the total number of grid points in each direction, and 
whether subgrid scale modelling is used or not, are given in 
Table [I] Each of these four models develops a first overtone 
pulsation with a period P of about 3.85 days (this is slightly 
shorter than one would expect if the structure of the entire 
star were contained in the computational domain). 



3.2 Time stepping 

If we consider the rate at which the cell average of the basic 
hydrodynamical variables p, I, and e changes as a function 
of time, the largest possible timestep is determined by the 
CFL-condition. But in practice all our simulations were lim- 
ited by the radiative timestep. 

In the interior regions the diffusion approximation is 
valid and thus 



Atdiff oc min 



/ 3 c p //tp 
V 16k<tT 3 VT 



(34) 



As usual, c p stands for the heat capacity at constant pressure 
and a for the Stefan-Boltzmann constant, k is essentially the 
inverse of the relevant, in our case therefore minimal, grid 
spacing, which enforces the most stringent time step restric- 
tion, i.e. k — C/ min(Ar^, Ayi) with typically 1 ^ C ^ 2n. 
C depends on the method used for spatial discretisation. In 
the code we use C = 1 and a constant safety factor less than 
or equal to 1 is applied in front of the entire expression. The 
resulting time scale is very small in the optically thin regions. 
On the other hand, the time scale for relaxing a temperature 
perturbation of arbitrary optical thickness by radiation (as- 
suming the radiative transfer equation where appropriate, 



not just the diffusion approximation) was shown by Spiegel 



( 1957 1 to be, in linear approximation of the perturbation, 



At, 



ad oc min 



16«:ctT 3 



KO Kp 

1 arccot — 

k k 



(35) 



Using radiative transfer in optically thin regions (and there- 
fore applying the timestep restrictions expressed by Eq. ( 35 1 



rather than Eq. ( 34 1 ) results in considerably larger allowed 
numerical timesteps for our Cepheid. Eq. ( 35 1 leads to a 
minimum for At ra( j in the optically thick regions where, how- 
ever, At ra d converges to the time scale of radiative diffusion 
anyway. The timestep restrictions for our explicit methods 
become ever more stringent when considering finer grids. 
Circumventing them requires an implicit method for the 
time integration of the Q ra d term which we have not yet 
implemented. 

3.3 Resolving steep gradients 

To resolve the steep gradients of density, temperature etc. 
in the hydrogen ionisation zone a resolution finer than 510 
points along the radial direction is needed. For a first es- 
timate of the required spacing we compare the pressure 
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We have checked that this estimate holds at least in ID and 
in Fig. [4] we present a comparison after a time evolution of 
0.25 sec from the initial, static model. Obviously in ID, a 
grid refinement of a factor of 2 is not enough, with a factor 
of 4, the radiative heating term starts to converge toward its 
actual shape and a factor of 8 suffices. The difference to a 
factor of 16 presently is not worth the extra computational 
costs required by such a fine resolution although it may still 
be large enough to lead to measurable effects in calculations 
of stellar spectra in future work. For our 2D simulations, 
where easing occurs due to the turbulent effects mentioned 
above, we applied a grid refinement factor of 3 in radial di- 
rection, and additional subgrid scale modelling. 



3 times 4times —2 times —coarse 



Figure 4. The shape of Q ra d (given in erg s 1 cm 3 ) for a 
one-dimensional simulation at different grid refinement factors. 



scale heights H p = p/|Vp| in the form p/(pg) for models of 
Cepheids and the Sun in the region of interest, i.e. the stel- 
lar photosphere at a Rosseland optical depth of order unity. 
Substituting the pressure of a perfect gas, p gas = pR gaa T//x, 
for the pressure p and the effective temperature T e ff for the 
temperature T thus yields a ratio of pressure scale heights 



T e ft M© 9e 



(36) 



Hp T e s M 9 
In our simulations of solar surface convection it was 
found that a minimum spacing of 20 km in radial direc- 
tion is necessary to just barely resolve all basic thermody- 
namical variables as a function of depth. This translates 
to approximately 4.25 Mm for our particular Cepheid. But 
we must also take into account the superadiabatic gradient 
V — Vad, because we have to resolve the temperature gradi- 
ent equally well as the pressure gradient. This may not be 



the case when considering only Eq. ( 36 1 to determine the 



spatial resolution of the simulation box, if the temperature 
gradient increases much more rapidly than it does at the su- 
peradiabatic peak of the solar convection zone. In the Sun 

Rosenthal et al.|1999[ >. For the 



max((V — V a 



0.6 (cf. 



initial condition of the Cepheid model, which was derived 
using a mixing-length treatment (MLT) of convection, the 
maximum temperature gradient critically depends on the 
a-parameter. In starting models where no turbulent pres- 
sure is included, max(V — V a d) ~ 14, thus the grid spacing 
would reduce to 0.18 Mm or a grid roughly seven times as 
fine as the grid used in model 1 at that location. On the 
other hand, if the starting model is derived with turbulent 
pressure included, max(V — V a d) ~ 3.2, and the grid of this 
model would only need to be refined roughly by a factor of 
two. While the first scenario is the only possible solution for 
a one dimensional, purely radiative model, there is reason 
to assume that in more dimensions a smaller factor can be 
achieved after convection has set in: the turbulent pressure 
caused by the flow supports the gas pressure in counteract- 
ing gravity and thus allows for lower gradients in gas pres- 
sure and temperature. The additional computation time is 
of course substantial: a factor of 7 in grid refinement leads 
to a factor ~ 20 for the number of temporal steps required 
when the timestep is limited by Af ra d as given by Eq. (351. 



4 GRID REFINEMENT 

For any study of the upper convection zone grid refine- 
ment is essential, as already pointed out in Sect. |3.3| To 
exclude effects that might result from fringing due to the 
grid refinement zone in an experiment with one dimensional 
simulations the gridspacing was reduced over the whole do- 
main. This improved the accuracy of Q ra d and, consequently, 
the lightcurve smoothed out. This can be seen in Fig. [5] 
which also shows the results from a comparison run with 
the lower, original resolution. Both simulations were con- 
tinued from the same initial model, an earlier state of the 
one-dimensional simulation on the coarser grid. Because of 
the higher spatial resolution the outgoing intensity becomes 
smaller: at the beginning the values are 30% above the 
average which is reached after an initial relaxation phase. 
The jagged lightcurve obtained with the coarse grid is to- 
tally unsuitable for comparisons with observations. This 
deficiency cannot be compensated by smoothing out the 
lightcurve since even then the average would remain too 
large in comparison with the simulation with lower resolu- 
tion and the shape of the curve would still be different from 
a sufficiently resolved calculation (cf. Fig. [5]l . The compu- 
tational resources required for the sufficiently resolved case 
run are unacceptably high for applications in two dimen- 
sions. Hence, the grid refinement has been modified to make 
it partially adaptive. The region of interest is the region 
where we expect to find the H-ionisation zone. Below this 
zone the maximum possible timestep decreases swiftly (see 
Fig. [6] for the inverse of Ai ra d), therefore a finer grid can 
be justified only above this border. Since the location of 
the H-ionisation zone varies both in angular direction and 
with time a fixed lower boundary would invariably lead to 
small timesteps, and thus we developed a co-moving lower 
boundary. The grid refinement is done over the entire angu- 
lar direction, since otherwise artifacts occur. 

The most important region for grid refinement is lo- 
cated around the maximum of the temperature gradient. To 
prevent the code from finding an inappropriate maximum 
we also make sure that the local temperature is at least T e g 
and the optical depth at least 100 at the lower border (black 
line in Fig. [Tf . 

We now briefly describe the refinement procedure. We 
consider a region where the maximum temperature gradi- 
ent can be found and the grid refinement factors gf (r) and 
gi(tp) are specified. If the maximum temperature gradient is 
outside this region, the region has to be reset. In the first 
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Figure 5. The light curve for a ID simulation on the coarse grid and on a grid having a 4 times higher resolution. Abscissa: output 
steps, ordinate: deviation from the average value. In both panels, the red line denotes the lightcurve while the black line indicates the 
radius. 




Figure 7. Increasing optical depth in the H ionisation zone on the fine grid (model 3). Approx. f° wide. White line: t : 
r = 10; black line t = 100; plotted over the convective flux F c = (u r — u r ) ■ I ph — ph) where h denotes the enthalpy. 



1; grey line: 



step the values of the physical quantities from the coarse 
grid are interpolated to the fine grid to obtain the start- 
ing state for the simulation. The step for the evolution in 
time is calculated for both grids. This yields the number 
N of steps on the fine grid during one step on the coarse 
grid. Because the calculations are limited by Ar. ra d, we get 
gf(r) ^ N ^ gf(r) 2 , where gf(r) is the radial grid refinement 
factor. Since the domains are not identical, it is possible that 
N is even smaller than gf(r), in which case we set N = gf(r). 
Now the "inner region" is determined. It reaches from the 
top of the maximum region to a few points below the lo- 
cation of the maximum of the temperature gradient. Grid 
refinement is only done in this "inner region". For the re- 
maining points near the superadiabatic peak the values are 
interpolated from the coarse grid. The lower boundary of 
the "inner region" is determined in every step, so that the 
bottom-line moves with the lower boundary of the actual H- 
ionisation zone. Artificial structures along this bottom line 
have not been observed by us in any of our simulations. 
Choosing, however, a top line in a similar manner we could 



not avoid getting artifacts: in the simulation shown in Fig. [8] 
such a moving top boundary was still used while this was 
no longer the case for the simulation shown in Fig. [9] Note 
that in all figures all values outside the "inner region" are 
only interpolated from the coarse grid. 

The time steps on the fine grid are done in exactly the 
same way as on the coarse grid. At the beginning of each 
step a time interpolation yields the values in the ghost cells 
at the top and bottom (black region in Fig. [8] and Fig. [9| 
of the grid refinement zone. In angular direction they are 
obtained from periodic continuation. After the last step in 
the grid refinement region the data are projected back to 
the coarse grid in a conservative way. 

For our two dimensional simulation (model 3) 3x4 ad- 
ditional grid points on each cell were used in addition to 
subgrid scale modelling. The resulting resolution in the crit- 
ical region is 0.66 Mm in radial and 1.16 Mm in angular 
direction leading to an aspect ratio of 1 : 1.8. But compar- 
ing to Fig. [4] admittedly there are still only 2 to 3 points to 
resolve the temperature gradient. 
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Figure 8. Grid refinement with upper limit on the "inner region". The left panel shows pressure, the right panel indicates the location 
of the grid refinement zone (grey: grid refinement zone, black: ghost cells). Note the structures in pressure (given in dyn cm -2 ) displayed 
in the left panel as indicated by arrows which are just at the location of the upper boundary indicated in the right panel. 




Figure 9. Grid refinement similar to Fig. [8] but this time without upper limit on the refinement zone ("inner region"). Artifacts no 
longer occur. 



The uppermost panel in Fig. [To] is taken from the pro- 
jection to the coarse grid and depicts the same area as the 
one below: the simulations both have an age of one day of 
stellar time. In both panels one can see the same maximum 
(red colour). However, not only is the convective flux already 
greater by a factor of 5, but there are also more convection 
cells in model 3 in comparison with the lower resolution 
model 2 and the convection cells are of an unnatural shape 
in model 2. Also, the maximum Mach numbers are signifi- 
cantly lower in model 2. The properties of the atmosphere 
and the differences brought about by resolution will be dis- 
cussed more closely in a subsequent paper. Already from 
these results, however, we can state that model 2 fails by far 
to properly represent the hydrogen ionisation zone and the 
higher resolution model 3 is needed. 



5 DISCUSSION 

5.1 Comparisons and applications of the models 

Among the simulations listed in Table [l] model 1 is partic- 
ularly useful for astrophysical studies. The simulation do- 
main is sufficiently wide such that the He n zone usually 
contains several up- and downflow structures. Its moderate 
resolution is both its main strength and shortcoming. Be- 
cause model 1 cannot resolve the hydrogen ionisation zone, 
it is very robust. Strong shock fronts develop much more 
rarely and, consequently, it has been easy to continue the 
simulation run over many pulsation cycles to provide sound 
statistical averages. Naturally, with this model studies on 
the convection-pulsation interaction have to be restricted to 
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Figure 10. Convective flux (in erg s — 1 cm -2 ) in the H-ionisation zone. The upper and lower panel depict the flux in the same area 
and at the same time according to the models with grid refinement (model 3, upper panel) and without (model 2, lower panel). The 
insufficiently resolved model 2 fails to properly develop the convection zone. Note that the same colours refer to drastically different 
values of the convective flux in both panels. Due to projection to the coarse grid for obtaining the upper panel, the refined model is not 
displayed in full resolution. 




Figure 6. Inverse of Ar ra( j, computed from the initial static 
model. Abscissa: gridpoints numbered from top to bottom, ordi- 
nate: C CO ur / A7"r ad in Hz. The H-ionisation zone is located at ap- 
proximately 100, the He II— ionisation zone at approximately 300. 



the role of the region around the zone of He II ionisation. 
An in-depth study based on this model which investigates 
enthalpy fluxes and the pressure work done as a function of 
the pulsation phase will be presented in paper II. 

Model 2 mainly serves as a reference simulation for 
model 3, which is derived from it through grid refinement. It 
has the same resolution as model 1 and differs from it only 
by a more narrow opening angle, thus it has no advantage 
over models 1 or model 3. 

On the other hand, model 3 and model 4 provide two 
alternative ways of performing high resolution simulations 
which do resolve the H I and He I ionisation zone while at the 
same time they include the He II ionisation zone and a suf- 
ficiently large region underneath it such that a meaningful 



and sufficiently accurate modelling of the radial pulsation is 
possible. As a consequence of its grid refinement model 3 has 
a wider range of applications than model 4: it offers about 
the same resolution combined with a much wider opening 
angle. Moreover, the maximum (resp. minimum) aspect ra- 
tio of the grid cells throughout the simulation domain can 
be reduced (kept closer to 1:1) while time step restrictions 
due to radiative diffusion in the layers below the superadia- 
batic layer can be eased as well. Finally, the He n ionisation 
zone in model 4 is affected by periodicity effects along the 
azimuthal coordinate which are caused by the necessarily 
small opening angle. This clearly demonstrates the advan- 
tages of grid refinement in combination with the spherical, 
radially stretched and co-moving mesh used in all our sim- 
ulations of Cepheids. From Fig. [10] it is evident that the 
dynamics at the optical surface of the Cepheid changes con- 
siderably once sufficient resolution is achieved: sharp, pro- 
nounced shock fronts appear which penetrate all the way 
up to the highest level of the photosphere. These phenom- 
ena will be discussed in detail in paper II. They appear in 
both model 3 and model 4 and the similarity of these fea- 
tures demonstrates that the grid refinement does not in- 
troduce any artifacts (the large, cusp like features at the 
optical surfaces which can be found in Fig. [8} {T0] are present 
in both model 3 and model 4 while they are hardly visible in 



model 2 in the lower panel of Fig. 10 1. Another consequence 



of successfully resolving the radiative cooling processes at 
the stellar surface is a much more realistic, relatively smooth 
lightcurve as shown by Fig. [5] for the one-dimensional case. 

However, the new level of dynamics found in high reso- 
lution models also introduces new numerical stability prob- 
lems. Indeed, for a long period of time we have been quite 
content with the closed boundary conditions at the top of 
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the domain as described in Sect. [277] That changed once the 
much more vigorous convection pattern in the H I and He I 
ionisation zone was found in the simulations with grid re- 



finement (see Fig. 10 1. When using the closed upper bound- 



ary condition as specified in Eq. ( 32 1 the faster flows and 



stronger shock waves generated in the highly resolved mod- 
els no longer lead to a benign behaviour as observed for the 
weaker shocks that appear in low resolution models. As a 
result, high resolution models cannot successfully be contin- 
ued over multiples of pulsation cycles. For such studies it is 
therefore necessary to develop appropriate variants of open 
boundary conditions. 



5.2 Overcoming the problem of small time steps 

To avoid unaffordable computational requirements and in- 
efficient use of computational resources high resolution sim- 
ulations of the hydrogen ionisation zone, as discussed with 
the examples of Fig. [4] for the one-dimensional case, require 
at least some terms to be treated with an implicit time inte- 
gration method. The feasibility of such an approach for the 
case of a moving or even an adaptive grid has been shown in 
Dorfi|(|1999|) for t he ID case. More recently, fViallet, Baraffe!] 
& Walder (20111 have developed an implicit 2D code and 



provided test cases for problems in stellar physics. 

These codes treat both the Euler equation and the dif- 
fusive term in the energy equation implicitly. For Cepheid 
modelling in multidimensions the situation is different from 
what the latter authors seem to envisage. What matters 
here is the stringent limitation which is imposed by the dif- 
fusive term in the energy on the timestep, if explicit time 
integration is used, rather than, e.g, restrictions due to the 
high sound speed near the bottom of the simulation domain. 
These limitations, which are caused by radiative diffusion, 
are due to unfavourable values of the physical parameters 
below the H+He I ionisation zone. On the other hand, the 
violent state of the atmosphere with strong shocks ensuing 
once the resolution is sufficient precludes, in this case, useful 
implicit time marching for the Euler equations themselves. 
That would be advantageous only in the case of fast flows 
with slowly moving features or for a low Mach number flow. 
For the numerical simulation of Cepheids it should hence 
be sufficient to implicitly advance the diffusive terms in the 
energy equation in time in the relevant regions. A necessary 
task in this respect is the development of novel Runge-Kutta 
methods which allow considerably larger time steps for such 
cases and replace the traditional schemes. Such new meth- 
ods have already been constructed and tested for their effi- 



ciency by Higueras et al. (20121, where their applicability is 



demonstrated for the case of semiconvection, although the 
new schemes have actually been developed with Cepheids 
and A-type stars in mind. 



6 CONCLUSIONS 

Setting the issue of strong shocks and obvious requirements 
such as the need for a moving grid and spherical geome- 
try aside Cepheid modelling in 2D and 3D is made difficult 
due to extremely different spatial scales which must be re- 
solved simultaneously. This problem becomes evident only 
once realistic microphysics, as characterised by the equation 



of state and radiative conductivities, combined with realistic 
values for fluxes and a realistic stratification due to gravity 
are taken into account. 

In our modified version of ANTARES we have tackled 
this challenge by introducing a grid in a spherical coordinate 
system stretched along the radial direction. The resolution is 
maximised by letting the grid move with the mean velocity 
of the uppermost domain layer and by introducing the possi- 
bility of grid refinement. The latter is essential to resolve the 
superadiabatic layer and the dynamical processes occurring 
in the photosphere such as sharp, fast moving shock fronts. 

These changes in turn necessitated modifications to the 
WENO scheme used for the spatial discretisation of the ad- 
vection operator and the pressure gradients in the dynamical 
equations. The new scheme discussed in this paper ensures 
conservation properties at the expense of a higher approx- 
imation order in smooth flow regions while it outperforms 
standard second order schemes due its smaller error con- 
stant. For the same reasons the discretisation of diffusion 
operators had to be modified accordingly. These modifica- 
tions allow to achieve a high spatial resolution at a fixed 
level of computational costs. 

A mandatory feature of a simulation code used to per- 
form realistic numerical simulations of Cepheids is a well 
scaling parallelisation. For this purpose the radiative trans- 
fer scheme used in ANTARES was not only adapted to the 
new grid, but also implemented in a new way such that com- 
putational nodes can be distributed efficiently both along ra- 
dial and azimuthal directions. This minimises overheads due 
to data communication and allows a cost-effective use of 256 
processor cores even with grids consisting only of typically 
around 500 by 500 points. 

Thanks to these features it has been possible to cal- 
culate a set of two-dimensional models which can be used 
to study the convection-pulsation interaction of Cepheids 
for realistic stellar parameters and realistic microphysics. A 
very high spatial resolution is needed to compute reliable 
light-curves and reliably calculate the stratification in the 
superadiabatic region. The grid refinement technique allows 
such resolution while covering a wide simulation domain and 
at the same time avoid the introduction of numerical arti- 
facts. 

For the time being two important problems remain on 
the technical and numerical side. Since unfavourably small 
timesteps are enforced by radiative diffusion below the H I 
and He I ionisation zone, the simulations are currently very 
expensive. A suitable implicit time integration procedure 
which eases these constraints is in development. Secondly, 
for high resolution simulations, open upper boundary con- 
ditions are required to allow stable integrations over many 
pulsation cycles. 

The models calculated thus far are already suited to in- 
vestigate the convection-pulsation interaction for the He II- 
ionisation zone with realistic microphysics and to get first in- 
sights to the atmospheric structure of Cepheids with (grey) 
radiative transfer from 2D models. Such investigations will 
be presented in paper II. 
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